Numerical Investigation of Glassy Dynamics in Low Density Systems 
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Vitrification in colloidal systems typically occurs at high densities driven by sharply varying, 
short-ranged interactions. The possibility of glassy behavior arising from smoothly varying, long- 
ranged particle interactions has received relatively little attention. Here we investigate the behavior 
of screened charged particles, and explicitly demonstrate that these systems exhibit glassy properties 
in the regime of low temperature and low density. Properties close to this low density (Wigner) 
glass transition share many features with their hard-sphere counterparts, but differ in quantitative 
aspects that may be accounted for via microscopic theoretical considerations. 

PACS numbers: 64.70.Pf, 61.20.Lc, 82.70.Dd 

The origins of the precipitous slowing down of dynam- 
ics in supercooled liquids are still unclear even after many 
decades of intense scrutiny. Model systems often form the 
basis for detailed investigations that include the core fea- 
tures known to give rise to generic glassy behavior. Per- 
haps the most prominent example of such a model system 
is the hard-sphere suspension, which has served as the 
basis for numerous experimental, theoretical and com- 
putational studies of the glass transition [l|, 0, 0|- Many 
of the most interesting properties of supercooled liquids 
and glasses, including two-step relaxation, stretched ex- 
ponential decay of density correlations and dynamic het- 
erogeneities occur in hard-sphere systems [1, 

While the glass transition of hard-spheres has become 
a paradigm for vitrification at high densities, serving as 
a reference point for conceptual attempts to connect the 
behavior of physically diverse classes of disordered ar- 
rested states of matter, another physically important 
limit of glassy systems has received much less system- 
atic scrutiny. This limit is that of a dilute assembly of 
particles interacting via long-ranged, soft repulsive forces 
such as those arising in charged systems. Over twenty 
years ago Chaikin and coworkers 0, Q investigated the 
phase behavior of dilute suspensions of charged colloids. 
Low density glasses stabilized by Coulomb repulsion were 
called "Wigner glasses" [ll- At that time a detailed inves- 
tigation of the structure and dynamics of these suspen- 
sions was not carried out. Hints of glassy behavior in 
one component plasmas have also been noted and the 
notion of a Wigner glass has been revived in colloidal 
systems 13, U, 12, 13] due to recent activity focusing on 
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gous to classical Wigner glass formation in colloidal sys- 
tems. 

Similar to the experimental situation, few theoretical 
investigations of the emergence of glassy dynamics in 
low density Coulomb systems have been carried out [2^. 
The most detailed investigations have been performed by 
Bosse and Wilke[23l. [23|, who have used idealized mode- 
coupling theory (MCT) jisl . to predict the dynamical 
behavior of a low density charged system in a neutral- 
izing background as it approaches the putative Wigner 
glass transition. These authors predict that glassy be- 
havior in this dilute regime shares many properties with 
the high density hard-sphere system, but some unique 
behavior also emerges. Indeed, Wigner glasses are stable 
even in the extreme dilute limit where the static struc- 
ture factor shows no modulation due to molecular shell 
structure. Furthermore, when the electrostatic repulsion 
is complemented by an excluded volume interaction, the 
MCT glass line is predicted to show a reentrant behavior 
due to the (density dependent) competition between the 
hard-core and the soft long-ranged repulsion 2^ 27 1. 

There are two main goals for the work presented 
here. First, we show that it is possible to generate 
Wigner glasses in simulations of dilute binary long ranged 
(screened) Coulomb system with a judicious choice of 
the interaction potential parameters and of the studied 
state point. Indeed, no previous theoretical work has ad- 
dressed the stability of the Wigner glass with resgect to 
the competing facility of crystallization [28i [29l [sol 31 1. 
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physical gelation in charged systems [IJ, [15!, 
Indeed, the study of glassy properties of dilute Coulomb 
systems has consequences that reach beyond classical 
systems and may shed light on routes to the formation 
of glasses in electronic systems 1^ 13 > where glassy ef- 
fects might persist even in the limit of weak to vanishing 
quenched disorder [2l|. Such self-induced glassiness re- 
sults from electron-electron interactions, an effect analo- 



This first result has important implications for the mod- 
eling and in terp retation of dynamics in charged colloidal 
suspensionsjUlill. Our second goal is to systematically 
study the glassy behavior of dilute Coulomb systems and 
to test some aspects of the MCT predictions for dynamics 
as the Wigner glass state is approached. 

Our choice of interparticle potential is dictated by 
a reasonable compromise between simplicity and real- 
ism. In particular, our simulations employ the stan- 
dard Yukawa form. Due to the rapid crystallization of 
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one-component Yukawa systems, we have investigated 
50% — 50% binary mixtures of point particles that may 
serve as models for dilute Wigner glasses. It should be 
noted that we have not included an additional short- 
range hard-sphere term in the potential. The advantage 
of this (unrealistic) simplification is that it allows us to 
study the effects of soft, long-ranged interactions in iso- 
lation of hard-core contributions. A disadvantage of this 
choice is that the notion of volume fraction is ill-defined, 
hence we will discuss our results in terms of number den- 
sity n = N/L^ where L is the simulation box length and 
N is the number of particles. The issue of reentrance, 
crucially depending on the existence of disparate length 
scales, may only be investigated if the hard-core repulsion 
is included and will be the subject of future work. 
The interaction potential between species i and j is. 



The amplitude of the repulsion differs for the two species 
as All = 0.2e, A22 = 3.5e, and A12 = 0.837e, with e 
being the unit of energy, while the screening lengths are 
all identical and taken to be Cii = ? = 1- Lengths are 
measured in units of ^. Such mixture corresponds to a 
model of two dilute colloidal species with different sur- 
face charges [3^. We perform molecular dynamics simula- 
tions of point particles of unitary mass m for two system 
syzes, i.e. N = 10"^ and N — 10^, at a fixed number 
density n = 0.002984, while the temperature T (mea- 
sured in units of ks, which we set equal to 1) is varied. 
In the following we always use the reduced temperature 
T* = r/10-5. The time step is 0.5 in units of^/{me)^/^. 
Since the potential is long-ranged, we choose a cutoff dis- 
tance equal to 25^, which ensures that all repulsive con- 
tributions at sufficiently long distances are accounted for. 
For all studied state points with T* > 1.7 data has been 
collected in the NVE ensemble following a sufficiently 
long equilibration period, controlling that no aging phe- 
nomena are detectable. Simulations for T* < 1.7 show 
aging even after extremely long equilibrations, confirm- 
ing the approach to a non-ergodic state. 

In Fig[l]we show results for the partial radial distribu- 
tion functions gijir) and for the partial static structure 
factors Sij{q) in the supercooled regime. A progressive 
structuring, as well as a shift towards larger distances, 
of the peaks of all partial radial distribution functions 
is observed upon decreasing T. Standard signatures of 
supercooled liquid structure, such as a split in the sec- 
ond peak, are evident. Type-1 particles are significantly 
closer on average than type-2 ones, due to the weaker 
repulsion. No sign of crystallization is noted during the 
whole duration of the runs. Focusing on the partial struc- 
ture factors, we notice that the amplitude of all Sij is 
small as compared to the canonical glass formers. More- 
over, Sii{q) shows a consistent upturn at low q, again a 
manifestation of the smaller repulsion experienced by the 





FIG. 1: Evolution in the supercooled regime, upon decreasing 
T, of: (a) partial radial distribution functions gii(r), (712 (r) 
and g22{r). The perfect agreement between two different sim- 
ulation sizes, i.e. N = 10^ (lines) and = lO" (symbols), 
provides evidence that results are equilibrated, reproducible 
and not phase-separating; (b) partial static structure factors 
Sij{q), as well as Sec (g) [331 which rules out the presence of a 
demixing transition. Symbols are results for simulations with 
— 10^ particles, lines are splined curves. 

type-1 particles with respect to the type-2 particles. In- 
deed, such feature is less evident for 5*12 while it is absent 
for 5*22 ■ To rule out the possibility that a phase separa- 
tion process interferes with the dynamical slowing down, 
we have monitored the time dependence of Sii{q) and 
visually inspected the configurations, both observations 
providing evidence against this putative process. More- 
over, upon further decrease of T, the low-g amplitude of 
Sii{q) saturates to a constant value. The concentration- 
concentration structure factor 5'cc('z)[3^, (shown in the 
inset of FiglUb)) does not display any increase at low q, 
ruling out the possibility of a demixing transition. Re- 
sults for the structural properties gij(r) and Sij{q) for 
the larger system size, also shown in Fig. [1] are identical 
to those of the smaller system. 

Having demonstrated that our system shows stable 
structural behavior, we next turn to investigate the dy- 
namical behavior. In Fig. [2]-a we show the decay of the 
self autocorrelation function of the density of type-1 par- 
ticles, 0f(g,<) close to the first peak of the structure 
factor. Strikingly, the behavior observed is qualitatively 
similar to that seen in the familiar dense systems contain- 
ing a harsh repulsive core. In particular, the dynamics 
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FIG. 2: Self autocorrelation function for the density of type- 
1 particles (j)i{q,t) at ~ 1.09: (a) T dependence (from left 
to right: T* = 10,5,3,2.5,2.25,2,1.9,1.7,1.5,1.3) and (b) 
TTS scaling only in the temperature range 1.7 ^ T* < 3. 
The dashed line is a stretched exponential fit of the master 
curve with stretching exponent f} 0.65. In (a) data for 
the larger simulation size (symbols) are also reported. Below 
T* = 1.7 the system shows a clear aging behavior, as well as a 
breakdown of TTS due to the increase in the plateau height. 
For type-2 particles, a similar, albeit slower dynamics (due to 
the larger repulsion amplitude), occurs. 



slows down by several orders of magnitude, giving rise to 
a typical two-step decay, in which the final relaxation is 
of stretched exponential form. Hence, the formation of 
a wide, soft cage around each particle arises, due to the 
long-range repulsive interactions of different strengths. 

In Fig. O-b we test time-temperature superposition 
(TTS) via rescahng cf)\{q,t) by the (T-dependent) a- 
relaxation time. For all 4'\{q,t) showing a two-step de- 
cay TTS holds down to T* = 1.7. For T* < 1.7 the 
plateau height changes significantly, leading to a clear 
breakdown of TTS. In the same T-range, aging effects 
can be observed, namely extremely slow drifts of the en- 
ergy as well as dependence of the dynamical variables of 
the initial observation times. 




FIG. 3: {r\{t)) for type-1 particles with decreasing T* . Lines 
and symbols refer respectively to simulations with A'' = 10"^ 
and A'^ = 10* particles. Inset: Power-law fit to the diffusion 
coefficient D, with exponent 7_d = 1.45 and T* ~ 1.67. 



We also study the mean squared displacement {rf{t)) 
as a function of T*. It shows the typical development 
of a long intermediate plateau, indicative of a dynami- 
cal slowing down (see Fig. The magnitude of the 
plateau is quite high, larger than ^. The extracted dif- 
fusion coefficient D, evaluated from the long time limit 
of the MSD, is reported in the inset. Power-law fits to 
D \T* — T*p° provide a method to estimate an effec- 
tive "mode-coupling temperature" of T* = 1.7, although 
the extracted exponent 7_d — 1-45 is smaller than the 
lowest limit set by MCT[25]. Below T*, we observe clear 
deviations from the power-law behavior suggesting that 
hopping processes are particularly relevant due the soft- 
ness of the repulsive cages. A further evidence of a change 
in the dynamics around T* is the fact that the plateau 
height, which can be opcratively defined as the inflection 
point oilog{MSD) vs. log(t), significantly decreases for 
T* < 1.7. 
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FIG. 4: Debye- Waller factors from simulations extracted 
from stretched exponential fits (symbols) and from MCT cal- 
culations (lines) for type-1 (left) and type-2 (right) particles 
respectively. The simulation temperature is T* — 2.0, while 
for MCT calculations we report results at the critical temper- 
ature r^fOT = 2.81. 

An interesting aspect of the MCT predictions [l^ is 
that the two-step relaxation is maintained even as struc- 
tural correlations induced by Coulomb repulsion are di- 
minished due to dilution effects. In particular, the Debye- 
Waller factors, = limt^oo(t>'-' {q.t) / Sjj{q) (where 
(j)^^ {q, t) is the collective intermediate scattering function) 
for both species are expected to show weak oscillations as 
a function of q that become still weaker for lower density, 
until oscillations are not observable ^j. While we could 
not access this ultra- low density limit, the Debye- Waller 
factors extracted from our simulations displayed in Fig.[3] 
show the precursor of this behavior. Indeed, the simula- 
tion data show weaker oscillations in , as compared to 
that seen at T* for the standard Kob- Andersen Lennard- 
Jones (KA) mixture [H, [s^l- Two other features are 
noteworthy in the context of this comparison. First, the 
range of wave vectors over which caging is present is much 
smaller than in standard high-density systems such as the 
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KA mixture. The small q values are reflective of caging 
at long length scales induced by the soft long-ranged po- 
tential. Secondly, ^ 1 for (7 ^ 0. This behavior is in 
contrast to that seen for the majority species in the KA 
mixture [stI, and might result from the relatively large 
compressibility of the system on these length scales [13|. 
It should be noted that both features are found in other 
soft materials such as colloidal gels(38j|. 

We also calculated /| directly from MCT using the 
numerical structure factors from Fig. [TJ Critical MCT 
Debye- Waller factors are also reported in Fig. [4] and they 
are compared to the ones calculated from fitting the nu- 
merical correlators for T* > T* . The agreement is re- 
markable, quantitatively matching all of the notable fea- 
tures exposed via direct MD simulation. 

Data presented in the previous figures show that a 
MCT description of the dynamics can be applied for 
T* > T* and that the qualitative relaxation behavior 
is of the usual "type-B" variety 25|, in agreement with 
the predictions of Bosse and Wilke 2^ 24 1 . While we 
postpone a detailed exposition for a future publication, 
it should be noted that features such the violation of the 
Stokes-Einstein relation, the growth of the non-Gaussian 
parameter and multi-point dynamical susceptibilities also 
occur in our simulated system. These features become 
prominent close to T* and behave in a manner similar to 
that found in standard glass-forming systems. 

In conclusion, we have clearly demonstrated the exis- 
tence of a stable classical Wigner glass state in a model 
of low-density charged particles, which is independent of 
competing processes such as crystallization or phase sep- 
aration. The low density soft glass studied here exhibits 
significant differences from the high-density hard-sphere 
glass. Most dramatically, we recall a much larger local- 
ization length, which manifests itself both in the MSD 
plateau as well as in the width of the Debye- Waller fac- 
tors. Despite these differences, several important as- 
pects of the dynamics can be successfully described by 
MCT. On the other hand, deviations from hard-sphere 
behaviour and from MCT predictions, in particular a 
small value of the power-law exponent for the diffusion 
constant as well as clear indications of activated pro- 
cesses, facilitated by the softness of the potential, are 
also observed. The possibility of reentrant relaxation em- 
anating from an interplay of hard core and long-ranged 
repulsion has not been explored in this work, and will be 
the subject of future investigation. 
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